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Fourier's Law from Schrodinger Dynamics 
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We consider a class of one-dimensional chains of weakly coupled many level systems. We present 
a theory which predicts energy diffusion within these chains for almost all initial states, if some 
concrete conditions on their Hamiltonians are met. By numerically solving the time dependent 
Schrodinger equation, we verify this prediction. Close to equilibrium we analyze this behavior in 
terms of heat conduction and compute the respective coefficient directly from the theory. 
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Almost two hundred years ago Fourier conjectured that 
temperature (or as we know today: energy) tends to dif- 
fuse through solids once close enough to equilibrium. His 
results may be stated in the following form 
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c(T)^T(x,t) 
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where u[T(x, i)] is the energy density at point x and 
time t, T the temperature, c(T) = the specific heat 
capacity and k the thermal conductivity. Despite the 
ubiquitous occurrence of this phenomenon, its explana- 
tion on the basis of some reversible microscopic dynamics 
remains a serious problem [lj, |2| • 

One approach to this subject is the Peierls-Boltzmann 
theory 3]. To explain the emergence of energy dif- 
fusion through an isolator based on quantum theory, 
Peierls essentially proposed a modified Boltzmann equa- 
tion. He replaced classical particles by quantized quasi- 
particles such as phonons and assumed statistical transi- 
tion rates taken from Fermi's Golden Rule for the colli- 
sion term. This concept faces conceptual shortcomings; 
though quantized normal modes are treated as classical 
particles, i.e., as being always well localized in configura- 
tion as well as momentum space. Furthermore, in order 
to exploit Fermi's Golden Rule in a classical picture, the 
complete actual quantum state of a phonon mode is dis- 
carded and only the mean occupation number is kept 
which is then treated as classical number of particles. 
Due to the neglect of any phases, this is called the ran- 
dom phase approximation and, as Peierls himself points 
out, so far lacks concrete justification. 

Another concept adressing the occurrence of regular 
heat conduction is the Green-Kubo formula. Derived on 
the basis of linear response theory it has originally been 
formulated for electrical transport 0, 0| . In this context 
the current is viewed as the response to a perturbative 
electrical potential which can be expressed as a part of 
the system Hamiltonian. But eventually the Green-Kubo 
formula boils down to a current-current auto-correlation, 
thus it can ad hoc be transfered to heat transport simply 
by replacing the electrical current by a heat current [fj. 
However, the justification of this replacement remains an 
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FIG. 1: Heat conduction model: N coupled subunits with 
ground level and band of n equally distributed levels (AE — 
1). Black dots refer to used initial states. 



open problem since there is no way of expressing a tem- 
perature gradient in terms of an addend to the Hamil- 
tonian [remarkably enough, Kubo himself comments on 
that replacement in a rather critical way Despite 
these unsolved questions, it has become a widely em- 
ployed technique 0, H, ■ 

To overcome this problem such Kubo-scenarios have 
recently been transfered from Hilbert- to Liouville space, 
where temperature gradients may be formulated in terms 
of operators [lOj . The method reveals normal heat trans- 
port |llHl2j | in very small quantum systems. However, it 
is numerically challenging especially for larger systems. 

In this letter we introduce yet another approach to heat 
transport within quantum systems, based on the Hilbert 
space Average Method (HAM) [II H [3. To demon- 
strate the emergence of heat diffusion directly from first 
principles (Schrodinger equation) we apply HAM to an 
appropriate "modular design model" , a weakly coupled 
chain of identical subsystems (cf. Fig.^l. HAM predicts 
that, if the model meets certain criteria (cf. (|14f> ) the lo- 
cal energy current between two adjacent subunits simply 
depends linearly on the difference of their inner energies. 
(For spatially small subunits: on the local gradient of the 
energy.) This diffusive behavior is, slightly reformulated 
[applying the continuity equation for the energy density 
S| = V • j to J[J] ! exactly what Fourier's law states 
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where j is the energy current density. To verify HAM we 
test its predictions for small systems (up to ten subunits) 
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FIG. 2: Probability to find the excitation in the fi = 1,2,3 
system. Comparison of the HAM prediction (lines) and the 
exact Schrodinger solution (dots). (N = 3, n = 500, A = 
5 ■ 10~ 5 , 5e = 0.05) 

by solving the corresponding time dependent Schrodinger 
equation numerically. HAM, however, allows for a dis- 
tinction between diffusive and non-diffusive behavior and 
produces the diffusion constant (cf. <|1U|) directly for ar- 
bitrarily large systems. 

The class of systems we are going to analyze is depicted 
in Fig. ^ with the Hamiltonian 

N N-l 

ff = ^ J ffioc(/i)+A^^M+l)- (3) 

f-i—i /i— i 

Here N identical subunits are assumed to have a non- 
degenerate ground state and a band of n exited states 
each, equally distributed over some band width Se such 
that the band width is small compared to the local energy 
gap AE (see Fig. [I] Se < AE, Se in units of AE). These 
subunits are coupled by an energy exchanging next neigh- 
bor interaction V((J,,fJ, + 1), chosen to be a (normalized) 
random hermitian matrix allowing for any possible tran- 
sition such as to avoid any bias. (Our results will turn 
out to be independent of the exact form of the matrix.) 
We choose the next neighbor coupling to be weak com- 
pared to the local gap (A <C AE, A in units of AE) . This 
way the full energy is approximately given by the sum of 
the local energies and those are approximately given by 
{iJj(t)\Hi oc (n)\ip(t)) = AEP^ where P p is the probabil- 
ity to find the /x'th subsystem in its excited state. This 
clean partition of the Hamiltonian into a strong local and 
a weak interaction part allows for a unique definition of 
both, a local energy (respectively temperature) as well as 
a local current. Again, this oversimplified model is pri- 
marily meant to demonstrate the possible direct emer- 
gence of diffusive behavior from Schrodinger dynamics. 
(For the impact of our results on more realistic systems 
see end of this letter.) 

The Hilbert space average method (HAM) is in essence 
a technique to produce a guess for the value of some quan- 
tity defined as a function of if itself is not known 
in full detail, only some features of it. Here it is used to 
produce a guess for some expectation value (ip\ A\ip) if the 
only information about |^>) is the given set of expectation 



values (tp\Pfj,\ip) — Pp. (Thus, P M is the projector pro- 
jecting out the subspace which corresponds to the /x'th 
subsystem occupying its excitation band, rather than its 
ground state.) Such a statement naturally has to be a 
guess, since there are in general many different \ip) that 
are in accord with the given set of P M but produce pos- 
sibly different values for (ip\A\ip) . The key question for 
the reliability of this guess, thus, is whether the distri- 
bution of the (V'|^4|V')' S produced by the respective set 
of |t/>)'s is broad or whether almost all those |^)'s yield 
(^>|.A|t/>)'s that are approximately equal. It can be shown 
that if the spectral width of A is not too large and A 
is high-dimensional almost all individual \if>) yield an ex- 
pectation value close to the mean of the distribution of 
the (^|4|V')' S an d thus the HAM guess will be reliable 
(see Fig. Eland furthermore jlij). To find that mean one 
has to average with respect to the \tl>)'s. This is called a 
Hilbert space average A and denoted as 

A = Tr{Aa} with a := M(1>\} {W p M)=Pll} ■ (4) 

This expression stands for the average of (ip\A\ip) — 
Tr{i|-0)(V>|} over all \tp) that feature {ip\P M \ip) = P M 
but are uniformly distributed otherwise. Uniformly 
distributed means invariant with respect to all uni- 
tary transformations U that leave P M unchanged, i.e., 
{ip\U^ PJj\i;) = P M . Thus the P's are specified by 
[U,Pfj] = 0. How is a to be computed? Any such U 
has to leave a invariant, i.e., [U, a] — 0. Furthermore, a 
has to obey TrjdP^} = P M . Those conditions uniquely 
determine a as 

-Ev^. (5) 

i.e., an expansion in terms of the projectors P^ them- 
selves. 

How can this be exploited to describe the dynamics 
of our system? Consider the full system's pure state at 
some time t, \tp(t)). Let D(t) be a time evolution op- 
erator describing the evolution of the system for a short 
time, i.e., \ip(t + r)) = D(T)\ip(t)) . This allows for the 
computation of the probabilities P M at time t + t finding 
the /x'th system somewhere in its excited band 

P M (* + r) = Tr{Z>(r) |V(t)>M*)| & (t) P,J . (6) 

Assume that rather than \ijj(t)) itself only the set of ex- 
pectation values P M is known. The application of HAM 
produces a guess for the P^(t + r) based on the P M (f) 
through replacing \tp(t))(tp(t)\ by the above a 

In the interaction picture the dynamics of the full system 
are only controlled by the interaction V{t). The time 
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evolution is generated by the corresponding Dyson series. 
Thus, assuming weak interactions Q may be evaluated 
to second order with respect to the interaction strength 
using an appropriately truncated Dyson series for D(t). 
This yields after extensive but rather straight forward 
calculations 



P^t + r) « fir) [P,-i{t) + P, +1 (t) - 2P M (t)] 



(8) 
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Tr{f(T")V r (0)}dr"dr' . (9) 



The above integrand is essentially the same environmen- 
tal correlation function that appears in the memory ker- 
nels of standard projection operator techniques. Those 
correlation functions typically feature some decay time 
r c after which they vanish. Thus, i.e, integrating them 
twice yields functions that increase linear in time after 
r c . Hence, for r > r c one simply gets /(r) = jt, where 7 
has to be computed from © but typically corresponds to 
a transition rate as obtained from Fermi's Golden Rule, 
i.e., 
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Assuming that the decay times of the correlation func- 
tions are small one can transform the iteration scheme 
(JSJ into a set of differential equations 



dPi 
dt 
dP M 
dt 
dP N 
dt 



-l(Pl ~ P2) , (11) 
- 7 (2P M - P M _x - P M+1 ) , (12) 

- 7 (Pv - Pv-i) • (13) 



Applying a discrete version of the continuity equation, 
dP^/dt = J M _i - J„, to (tTTlTSI) yields J M := -7(^+1 - 
P M ). This reduces to (J2J), i.e., Fourier's law, with 7 = k/c. 

The above scheme, however, only applies if the dynam- 
ics of the system are reasonably well described by a Dyson 
series truncated at second order also for times r larger 
than t c . This and similar arguments yield the following 
necessary conditions for the above described occurrence 
of diffusive transport 



(14) 



We checked by numerically integrating the Schrodinger 
equation (see below) that no diffusive transport as de- 
scribed by H11I13JI results if those criteria are violated. 
(For a more detailed description of HAM and its various 
implications, see 0,0].) 

To analyze validity and performance of HAM we com- 
pare its results with data from a direct numerical integra- 
tion of the Schrodinger equation. This is, of course, only 
possible for systems small enough to allow for the latter. 
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FIG. 3: Deviation of HAM from the exact solution {D\ least 
squares): dependence of D\ on n for N = 3. Inset: depen- 
dence of D\ on N for n = 500. 



Hereby we restrict ourselves to initial states with only one 
subsystem in the exited band (all others in the ground 
level, black dots in Fig.QJ. Finding an effective Hamilto- 
nian for the one-excitation subspace we are able to solve 
the Schrodinger equation for up to N = 10 subsystems, 
n = 500 levels each. Firstly restricting to N = 3, the nu- 
merical results together with the HAM predictions for an 
initial state with Pi = 1, Pj = P3 = are shown in Fig. [21 
(N = 3, n = 500, 5e = 0.05, A = 5 • 10~ 5 ). There is a rea- 
sonably good agreement. To investigate the accuracy of 
the HAM for, e.g., P\{t) we introduce Df, being the time- 
averaged quadratic deviation of the exact (Schrodinger) 
result for Pi (t) from the HAM prediction 



(P* AM (t) - Pl™ ct {t) fdt . (15) 



To analyze how big this "typical deviation" is, we have 
computed D\ for a N — 3 system with the first subunit 
initially in an arbitrary excited state, for different num- 
bers of states n in the bands. As shown in Fig. the 
deviation scales like 1/n with the band size, i.e., vanishes 
in the limit of high dimensional subunits. This behavior 
does not come as a complete surprise since it has the- 
oretically been conjectured and numerically verified in 
the context of equilibrium fluctuations of non-Markovian 
systems The inset of Fig. |3] shows that D\ goes 

down also with increasing number of subunits N but then 
levels off. Altogether HAM appears to be applicable even 
down to moderately sized systems. So far we have re- 
stricted ourselves to pure states. A drastic further re- 
duction of D 2 can be expected for mixed states (which 
are typical in the context of thermodynamical phenom- 
ena), since pure state fluctuations can be expected to 
cancel partially if added together. 

So far we have considered energy diffusion through the 
system only. The final state should approach equiparti- 
tion of energy over all subunits (see Fig. |5Jl - a thermal 
equilibrium state (see 0|). Close to this equilibrium we 
expect the system to be in a state where the probability 
distribution of each subunit is approximately canonical 
(Gibbs state) but still with a slightly different tempera- 
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FIG. 4: Heat conductivity (1161 over temperature for a system 
with n = 500, Se = 0.05 and A = 5 • 10" 5 . 



ture T M for each site (see Fig. |2 for t > 1000). Special- 
izing in those "local equilibrium states" and exploiting 
the HAM results allows for a direct connection of 

the local energy current between any two adjacent sub- 
units with their temperature gradient AT = T\ — T% and 
their mean temperature T = (T\ + T2)/2. Since this con- 
nection is found to be linear in the temperature gradient 
one can simply read off the temperature dependent heat 
conductivity (cf. J2J) 
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as displayed in Fig. 0] For our (simple) model this is in 
agreement with k = jc (as stated after pi!13|l ). if one 
inserts for c the specific heat of one subunit. The result 
of Fig. 0| is similar to the thermal conductivity of gaped 
quantum spin chains as calculated by the "Liouville space 
method" 

What impact do those results for our "design model" 
have on real physical systems? For an application of 
HAM one has to organize the system as a "net-structure" 
of weakly coupled subunits in the sense of ((HI) . If this 
can be established HAM predicts that energy will dif- 
fuse from subunit to subunit, irrespective of whether the 
coupled subunits form a one- or a more-dimensional net- 
structure. Such a structure may be achieved by coarse 
graining periodic systems like spin chains, crystals, etc. 
in entities containing many elementary cells and taking 
the mesoscopic entities as the subunits proper and the ef- 
fective interactions between the entities as the couplings. 
This way increasing the "grain size" will result in higher 
state densities within the subunits and relatively weaker 
couplings such that above minimum grain size the criteria 
(fHtjl may eventually be fulfilled. 

Of course, the resulting subunits cannot generally be 
expected to feature the same gaped spectral structure as 
our design model. (Although there are, e.g., spin sys- 
tems which do so 0,0]-) HAM also applies to multi- 



band subunits. Then, for normal transport, the bands 
and their couplings that mediate the transport should 
obey the criterion p4|) individually. If this is the case 
the gaps between the bands are dispensable. Only the 
band width has to be small compared to the mean band 
energy. Thus, even a continuous spectrum may be spec- 
trally coarse grained to fulfill (|14fl . This generally yields 
different energy diffusion constants 7 for different bands. 
Nevertheless a concrete temperature dependent conduc- 
tivity k would result for close to equilibrium states. Thus 
HAM should be applicable also to more realistic systems, 
and work in that direction is under way. 

In conclusion we have shown that already a simple 
modular quantum system can give rise to diffusive en- 
ergy transport. These results are attained based on a 
new efficient approximation scheme and confirmed by a 
full Schrodinger analysis. 
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